package hos.map.tool;

import java.util.List;

import hos.map.coordinate.LatLng;

/**
 * <p>Title: PolygonUtil </p>
 * <p>Description:  </p>
 * <p>Company: www.mapuni.com </p>
 *
 * @author : 蔡俊峰
 * @version : 1.0
 * @date : 2023-12-07 22:02
 */
public class PolygonUtil {

    private PolygonUtil() {
    }

    private static double getEarthRadius() {
        return MathUtil.getEarthRadius();
    }

    private static double getPI() {
        return MathUtil.getPI();
    }

    /// Returns tan(latitude-at-lng3) on the great circle (lat1, lng1) to
    /// (lat2, lng2). lng1==0.
    /// See http://williams.best.vwh.net/avform.htm .
    private static double _tanLatGC(double lat1, double lat2, double lng2, double lng3) {
        return (Math.tan(lat1) * Math.sin(lng2 - lng3) + Math.tan(lat2) * Math.sin(lng3)) / Math.sin(lng2);
    }


    /// Returns mercator(latitude-at-lng3) on the Rhumb line (lat1, lng1) to
    /// (lat2, lng2). lng1==0.
    private static double _mercatorLatRhumb(double lat1, double lat2, double lng2, double lng3) {
        return (MathUtil.mercator(lat1) * (lng2 - lng3) + MathUtil.mercator(lat2) * lng3) / lng2;
    }

    /**
     * 计算到South Pole的垂直线段(lat3, lng3)是否与到(lat1, lng1)的线段(lat2, lng2)相交。
     * 经度偏移量为- long1;隐式的ln1变成0。
     * /// Computes whether the vertical segment (lat3, lng3) to South Pole
     * /// intersects the segment (lat1, lng1) to (lat2, lng2).
     * /// Longitudes are offset by -lng1; the implicit lng1 becomes 0.
     */
    private static boolean _intersects(double lat1, double lat2, double lng2, double lat3, double lng3, boolean geodesic) {
        // Both ends on the same side of lng3.
        if ((lng3 >= 0 && lng3 >= lng2) || (lng3 < 0 && lng3 < lng2)) {
            return false;
        }
        double pi = getPI();
        // Point is South Pole.
        if (lat3 <= -pi / 2) {
            return false;
        }
        // Any segment end is a pole.
        if (lat1 <= -pi / 2 || lat2 <= -pi / 2 || lat1 >= pi / 2 || lat2 >= pi / 2) {
            return false;
        }
        if (lng2 <= -pi) {
            return false;
        }

        final double linearLat = (lat1 * (lng2 - lng3) + lat2 * lng3) / lng2;
        // Northern hemisphere and point under lat-lng line.
        if (lat1 >= 0 && lat2 >= 0 && lat3 < linearLat) {
            return false;
        }
        // Southern hemisphere and point above lat-lng line.
        if (lat1 <= 0 && lat2 <= 0 && lat3 >= linearLat) {
            return true;
        }
        // North Pole.
        if (lat3 >= pi / 2) {
            return true;
        }

        // Compare lat3 with latitude on the GC/Rhumb segment corresponding to lng3.
        // Compare through a strictly-increasing function (tan() or
        // MathUtil.mercator()) as convenient.
        return geodesic ? Math.tan(lat3) >= _tanLatGC(lat1, lat2, lng2, lng3) : MathUtil.mercator(lat3) >= _mercatorLatRhumb(lat1, lat2, lng2, lng3);
    }

    /**
     * 计算给定点是否在指定多边形内。
     */
    public static boolean containsLocation(LatLng point, List<LatLng> polygon, boolean geodesic) {
        return containsLocationAtLatLng(point.getLatitude(), point.getLongitude(), polygon, geodesic);
    }

    /**
     * 计算给定点是否在指定多边形内。多边形总是被认为是封闭的，不管最后一个点是否等于第一个点。
     * 内部被定义为不包含南极——南极总是在外部。如果测地线为真，则多边形由大圆段组成，否则由菱形(凸形)段组成。
     * /// Computes whether the given point lies inside the specified polygon.
     * /// The polygon is always considered closed, regardless of whether the last
     * /// point equals the first or not.
     * /// Inside is defined as not containing the South Pole -- the South Pole is
     * /// always outside. The polygon is formed of great circle segments if geodesic
     * /// is true, and of rhumb (loxodromic) segments otherwise.
     */
    public static boolean containsLocationAtLatLng(double latitude, double longitude, List<LatLng> polygon, boolean geodesic) {
        if (polygon == null || polygon.isEmpty()) {
            return false;
        }

        final double lat3 = MathUtil.toRadians(latitude);
        final double lng3 = MathUtil.toRadians(longitude);
        final LatLng prev = polygon.get(polygon.size() - 1);
        double lat1 = MathUtil.toRadians(prev.getLatitude());
        double lng1 = MathUtil.toRadians(prev.getLongitude());
        int nIntersect = 0;
        double pi = getPI();
        for (LatLng point2 : polygon) {
            final double dLng3 = MathUtil.wrap(lng3 - lng1, -pi, pi);
            // Special case: point equal to vertex is inside.
            if (lat3 == lat1 && dLng3 == 0) {
                return true;
            }
            final double lat2 = MathUtil.toRadians(point2.getLatitude());
            final double lng2 = MathUtil.toRadians(point2.getLongitude());
            // Offset longitudes by -lng1.
            if (_intersects(lat1, lat2, MathUtil.wrap(lng2 - lng1, -pi, pi), lat3, dLng3, geodesic)) {
                ++nIntersect;
            }
            lat1 = lat2;
            lng1 = lng2;
        }
        return (nIntersect & 1) != 0;
    }

    private static final double defaultTolerance = 0.1; // meters.

    /**
     * 在以米为单位的指定公差范围内，计算给定点是否位于多边形的边缘上或边缘附近。
     * 如果测地线为真，则多边形边缘由大圆段组成，否则由菱形段组成。多边形边缘是隐式封闭的-包括第一个点和最后一个点之间的封闭段。
     * /// Computes whether the given point lies on or near the edge of a polygon,
     * /// within a specified tolerance in meters. The polygon edge is composed of
     * /// great circle segments if geodesic is true, and of Rhumb segments
     * /// otherwise. The polygon edge is implicitly closed - the closing segment
     * /// between the first point and the last point is included.
     */
    public static boolean isLocationOnEdge(LatLng point, List<LatLng> polygon, boolean geodesic) {
        return isLocationOnEdge(point, polygon, geodesic, defaultTolerance);
    }

    /**
     * 在以米为单位的指定公差范围内，计算给定点是否位于多边形的边缘上或边缘附近。
     * 如果测地线为真，则多边形边缘由大圆段组成，否则由菱形段组成。多边形边缘是隐式封闭的-包括第一个点和最后一个点之间的封闭段。
     * /// Computes whether the given point lies on or near the edge of a polygon,
     * /// within a specified tolerance in meters. The polygon edge is composed of
     * /// great circle segments if geodesic is true, and of Rhumb segments
     * /// otherwise. The polygon edge is implicitly closed - the closing segment
     * /// between the first point and the last point is included.
     */
    public static boolean isLocationOnEdge(LatLng point, List<LatLng> polygon, boolean geodesic, double tolerance) {
        return _isLocationOnEdgeOrPath(point, polygon, true, geodesic, tolerance);
    }

    /**
     * 计算给定点是否位于折线上或附近，在指定的公差(以米为单位)内。如果测地线为真，
     * 则折线由大圆段组成，否则由菱形段组成。折线不闭合——第一个点和最后一个点之间的闭合段不包括在内。
     * /// Computes whether the given point lies on or near a polyline, within a
     * /// specified tolerance in meters. The polyline is composed of great circle
     * /// segments if geodesic is true, and of Rhumb segments otherwise.
     * /// The polyline is not closed -- the closing segment between the first point
     * /// and the last point is not included.
     */
    public static boolean isLocationOnPath(LatLng point, List<LatLng> polyline, boolean geodesic) {
        return isLocationOnPath(point, polyline, geodesic, defaultTolerance);
    }

    /**
     * 计算给定点是否位于折线上或附近，在指定的公差(以米为单位)内。如果测地线为真，
     * 则折线由大圆段组成，否则由菱形段组成。折线不闭合——第一个点和最后一个点之间的闭合段不包括在内。
     * /// Computes whether the given point lies on or near a polyline, within a
     * /// specified tolerance in meters. The polyline is composed of great circle
     * /// segments if geodesic is true, and of Rhumb segments otherwise.
     * /// The polyline is not closed -- the closing segment between the first point
     * /// and the last point is not included.
     */
    public static boolean isLocationOnPath(LatLng point, List<LatLng> polyline, boolean geodesic, double tolerance) {
        return _isLocationOnEdgeOrPath(point, polyline, false, geodesic, tolerance);
    }


    private static boolean _isLocationOnEdgeOrPath(LatLng point, List<LatLng> poly, boolean closed, boolean geodesic, double toleranceEarth) {
        final int idx = locationIndexOnEdgeOrPath(point, poly, closed, geodesic, toleranceEarth);

        return idx >= 0;
    }

    /**
     * 在指定的公差范围内计算给定点是否位于折线上或折线附近(以及在何处)。折线不闭合——第一个点和最后一个点之间的闭合段不包括在内。
     * /// Computes whether (and where) a given point lies on or near a polyline,
     * /// within a specified tolerance. The polyline is not closed -- the closing
     * /// segment between the first point and the last point is not included.
     * ///
     * /// @param point     our needle
     * /// @param poly      our haystack
     * /// @param geodesic  the polyline is composed of great circle segments if
     * ///                  geodesic is true, and of Rhumb segments otherwise
     * /// @param tolerance tolerance (in meters)
     * /// @return -1 if point does not lie on or near the polyline.
     * /// 0 if point is between poly[0] and poly[1] (inclusive),
     * /// 1 if between poly[1] and poly[2],
     * /// ...,
     * /// poly.size()-2 if between poly[poly.size() - 2] and poly[poly.size() - 1]
     */
    public static int locationIndexOnPath(LatLng point, List<LatLng> poly, boolean geodesic) {
        return locationIndexOnEdgeOrPath(point, poly, false, geodesic, defaultTolerance);
    }

    /**
     * 在指定的公差范围内计算给定点是否位于折线上或折线附近(以及在何处)。折线不闭合——第一个点和最后一个点之间的闭合段不包括在内。
     * /// Computes whether (and where) a given point lies on or near a polyline,
     * /// within a specified tolerance. The polyline is not closed -- the closing
     * /// segment between the first point and the last point is not included.
     * ///
     * /// @param point     our needle
     * /// @param poly      our haystack
     * /// @param geodesic  the polyline is composed of great circle segments if
     * ///                  geodesic is true, and of Rhumb segments otherwise
     * /// @param tolerance tolerance (in meters)
     * /// @return -1 if point does not lie on or near the polyline.
     * /// 0 if point is between poly[0] and poly[1] (inclusive),
     * /// 1 if between poly[1] and poly[2],
     * /// ...,
     * /// poly.size()-2 if between poly[poly.size() - 2] and poly[poly.size() - 1]
     */
    public static int locationIndexOnPath(LatLng point, List<LatLng> poly, boolean geodesic, double tolerance) {
        return locationIndexOnEdgeOrPath(point, poly, false, geodesic, tolerance);
    }

    /**
     * 在指定的公差范围内计算给定点是否位于折线上或折线附近(以及在何处)。如果关闭，则不考虑折线的最后一个点和第一个点之间的关闭段。
     * /// Computes whether (and where) a given point lies on or near a polyline,
     * /// within a specified tolerance. If closed, the closing segment between the
     * /// last and first points of the polyline is not considered.
     * ///
     * /// @param point          our needle
     * /// @param poly           our haystack
     * /// @param closed         whether the polyline should be considered closed by
     * ///                       a segment connecting the last point back to the
     * ///                       first one
     * /// @param geodesic       the polyline is composed of great circle segments if
     * ///                       geodesic is true, and of Rhumb segments otherwise
     * /// @param toleranceEarth tolerance (in meters)
     * /// @return -1 if point does not lie on or near the polyline.
     * /// 0 if point is between poly[0] and poly[1] (inclusive),
     * /// 1 if between poly[1] and poly[2],
     * /// ...,
     * /// poly.size()-2 if between poly[poly.size() - 2] and poly[poly.size() - 1]
     */
    public static int locationIndexOnEdgeOrPath(LatLng point, List<LatLng> poly, boolean closed, boolean geodesic, double toleranceEarth) {
        if (poly == null || poly.isEmpty()) {
            return -1;
        }
        double earthRadius = getEarthRadius();
        double pi = getPI();
        final double tolerance = toleranceEarth / earthRadius;
        final double havTolerance = MathUtil.hav(tolerance);
        final double lat3 = MathUtil.toRadians(point.getLatitude());
        final double lng3 = MathUtil.toRadians(point.getLongitude());
        final LatLng prev = closed ? poly.get(poly.size() - 1) : poly.get(0);
        double lat1 = MathUtil.toRadians(prev.getLatitude());
        double lng1 = MathUtil.toRadians(prev.getLongitude());
        int idx = 0;
        if (geodesic) {
            for (LatLng point2 : poly) {
                final double lat2 = MathUtil.toRadians(point2.getLatitude());
                final double lng2 = MathUtil.toRadians(point2.getLongitude());
                if (_isOnSegmentGC(lat1, lng1, lat2, lng2, lat3, lng3, havTolerance)) {
                    return Math.max(0, idx - 1);
                }
                lat1 = lat2;
                lng1 = lng2;
                idx++;
            }
        } else {
            // We project the points to mercator space, where the Rhumb segment is a
            // straight line, and compute the geodesic distance between point3 and the
            // closest point on the segment. This method is an approximation, because
            // it uses "closest" in mercator space which is not "closest" on the
            // sphere -- but the error is small because "tolerance" is small.
            final double minAcceptable = lat3 - tolerance;
            final double maxAcceptable = lat3 + tolerance;
            double y1 = MathUtil.mercator(lat1);
            final double y3 = MathUtil.mercator(lat3);
            final double[] xTry = new double[3];
            for (LatLng point2 : poly) {
                final double lat2 = MathUtil.toRadians(point2.getLatitude());
                final double y2 = MathUtil.mercator(lat2);
                final double lng2 = MathUtil.toRadians(point2.getLongitude());
                if (Math.max(lat1, lat2) >= minAcceptable && Math.min(lat1, lat2) <= maxAcceptable) {
                    // We offset longitudes by -lng1; the implicit x1 is 0.
                    final double x2 = MathUtil.wrap(lng2 - lng1, -pi, pi);
                    final double x3Base = MathUtil.wrap(lng3 - lng1, -pi, pi);
                    xTry[0] = x3Base;
                    // Also explore MathUtil.wrapping of x3Base around the world in both
                    // directions.
                    xTry[1] = x3Base + 2 * pi;
                    xTry[2] = x3Base - 2 * pi;
                    for (double x3 : xTry) {
                        final double dy = y2 - y1;
                        final double len2 = x2 * x2 + dy * dy;
                        final double t = len2 <= 0 ? 0 : MathUtil.clamp((x3 * x2 + (y3 - y1) * dy) / len2, 0, 1);
                        final double xClosest = t * x2;
                        final double yClosest = y1 + t * dy;
                        final double latClosest = MathUtil.inverseMercator(yClosest);
                        final double havDist = MathUtil.havDistance(lat3, latClosest, x3 - xClosest);
                        if (havDist < havTolerance) {
                            return Math.max(0, idx - 1);
                        }
                    }
                }
                lat1 = lat2;
                lng1 = lng2;
                y1 = y2;
                idx++;
            }
        }
        return -1;
    }

    /// Returns sin(initial bearing from (lat1,lng1) to (lat3,lng3) minus initial
    /// bearing from (lat1, lng1) to (lat2,lng2)).
    private static double _sinDeltaBearing(double lat1, double lng1, double lat2, double lng2, double lat3, double lng3) {
        final double sinLat1 = Math.sin(lat1);
        final double cosLat2 = Math.cos(lat2);
        final double cosLat3 = Math.cos(lat3);
        final double lat31 = lat3 - lat1;
        final double lng31 = lng3 - lng1;
        final double lat21 = lat2 - lat1;
        final double lng21 = lng2 - lng1;
        final double a = Math.sin(lng31) * cosLat3;
        final double c = Math.sin(lng21) * cosLat2;
        final double b = Math.sin(lat31) + 2 * sinLat1 * cosLat3 * MathUtil.hav(lng31);
        final double d = Math.sin(lat21) + 2 * sinLat1 * cosLat2 * MathUtil.hav(lng21);
        final double denom = (a * a + b * b) * (c * c + d * d);
        return denom <= 0 ? 1 : (a * d - b * c) / Math.sqrt(denom);
    }

    private static boolean _isOnSegmentGC(double lat1, double lng1, double lat2, double lng2, double lat3, double lng3, double havTolerance) {
        final double havDist13 = MathUtil.havDistance(lat1, lat3, lng1 - lng3);
        if (havDist13 <= havTolerance) {
            return true;
        }
        final double havDist23 = MathUtil.havDistance(lat2, lat3, lng2 - lng3);
        if (havDist23 <= havTolerance) {
            return true;
        }
        final double sinBearing = _sinDeltaBearing(lat1, lng1, lat2, lng2, lat3, lng3);
        final double sinDist13 = MathUtil.sinFromHav(havDist13);
        final double havCrossTrack = MathUtil.havFromSin(sinDist13 * sinBearing);
        if (havCrossTrack > havTolerance) {
            return false;
        }
        final double havDist12 = MathUtil.havDistance(lat1, lat2, lng1 - lng2);
        final double term = havDist12 + havCrossTrack * (1 - 2 * havDist12);
        if (havDist13 > term || havDist23 > term) {
            return false;
        }
        if (havDist12 < 0.74) {
            return true;
        }
        final double cosCrossTrack = 1 - 2 * havCrossTrack;
        final double havAlongTrack13 = (havDist13 - havCrossTrack) / cosCrossTrack;
        final double havAlongTrack23 = (havDist23 - havCrossTrack) / cosCrossTrack;
        final double sinSumAlongTrack = MathUtil.sinSumFromHav(havAlongTrack13, havAlongTrack23);
        return sinSumAlongTrack > 0; // Compare with half-circle == pi using sign of sin().
    }


    /**
     * 如果提供的点列表是封闭多边形(即第一个点和最后一个点是相同的)则返回true，如果不是则返回false
     * /// Returns true if the provided list of points is a closed polygon (i.e., the
     * /// first and last points are the same), and false if it is not
     * ///
     * /// @param poly polyline or polygon
     * /// @return true if the provided list of points is a closed polygon (i.e., the
     * /// first and last points are the same), and false if it is not
     */
    public static boolean isClosedPolygon(List<LatLng> poly) {
        return poly.get(0).equals(poly.get(poly.size() - 1));
    }

}
